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Abstract 

Monte Carlo (MC) techniques are often used to estimate integrals of a multivariate function using randomly gen- 
erated samples of the function. In light of the increasing interest in uncertainty quantification and robust design ap- 
plications in aerospace engineering, the calculation of expected values of such functions (e.g. performance measures) 
becomes important. However, MC techniques often suffer from high variance and slow convergence as the number 
of samples increases. In this paper we present Stacked Monte Carlo (StackMC), a new method for post-processing an 
existing set of MC samples to improve the associated integral estimate. StackMC is based on the supervised learning 
techniques of fitting functions and cross validation. It should reduce the variance of any type of Monte Carlo integral 
estimate (simple sampling, importance sampling, quasi-Monte Carlo, MCMC, etc.) without adding bias. We report 
on an extensive set of experiments confirming that the StackMC estimate of an integral is more accurate than both the 
associated unprocessed Monte Carlo estimate and an estimate based on a functional fit to the MC samples. These ex- 
periments run over a wide variety of integration spaces, numbers of sample points, dimensions, and fitting functions. 
In particular, we apply StackMC in estimating the expected value of the fuel burn metric of future commercial aircraft 
and in estimating sonic boom loudness measures. We compare the efficiency of StackMC with that of more standard 
methods and show that for negligible additional computational cost significant increases in accuracy are gained. 
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1 Introduction 

A system optimized only for best performance at nominal operating conditions can see severe performance degradation 
at even slightly off-design conditions. Aerospace systems rarely operate at precisely the design condition and a robust 
design approach dictates trading some performance at the nominal condition for improved performance over a wide 
range of input conditions. However, certain input conditions will occur rarely or never, and adding robustness for these 
conditions will degrade average system performance. Consequently, we are often interested in optimizing the expected 
performance of a system rather than the performance at any specific operating point. 
Formally, we are interested in an integral of the form 

B[/]=/= J f(x)p(x)dx , (1) 

where f(x) is the (potentially multivariate) objective function to be optimized, and p(x) is the probability density 
function (or probability distribution, in which case ([1} is actually a summation) from which x is generated. 

Evaluating ([TJ is non-trivial, especially when the objective function is high dimensional and/or expensive to eval- 
uate. When standard quadrature techniques (such as Simpson's rule) become intractable, Monte Carlo (MC) methods 
are powerful techniques which have become the standard for integral estimation. However, MC techniques also have 
their drawbacks; the required computational expense may be prohibitive in most situations. The question is if a tech- 
nique exists to significantly lower the cost of estimating (fTJ. This paper describes such a technique which combines 
MC with machine learning and statistical techniques and can lead to significant computational savings over traditional 
methods. 

We begin this paper with a brief review of integral estimation techniques including Monte Carlo, fitting algorithms, 
and stacking. We then introduce a new technique we call Stacked Monte Carlo (StackMC), which uses stacking 
to reduce the estimation error of any Monte Carlo method. Finally, we apply StackMC to a number of analytic 
example problems and to two problems from the aerospace literature. We show that StackMC can significantly reduce 
estimation error and never has higher estimation error than the best of Monte Carlo and the chosen fitting algorithm. 

2 Integral Estimation Techniques 
2.1 Estimation Error 

When using any method M returning /m as an estimate of Q, there are two sources of estimation error: bias (b) and 
variance (v). Bias is the expected difference between the method and the truth: b = E[/m - /], where the expectation 
is over all data sets. A method whose average over many runs gives the correct answer has zero bias, whereas a 



2 



method that estimates high or low on average has non-zero bias. As an example, the Euler equations have significant 
bias in their estimate of viscous drag. Variance is a measure of how much fu varies between different runs of M: 
v = E[(/m - E[/m]) 2 ]. If M is deterministic, it has zero variance because every run returns the same answer, whereas 
if M is stochastic multiple runs have different outputs leading to a non-zero variance. The total expected squared error 
is the combination of these two factors 

E[|/m-/T] =b 2 + v . 
Any method seeking to reduce estimation error must keep both of sources of error low. 

2.2 Monte Carlo Techniques 

In "Simple Sampling" Monte Carlo (MC), a set of N samples, D x , is generated randomly according to p(x). The 
estimate of the expected value of the function based on D x is the average of the function values of the samples 

1 N 

f*fmc = ^Y J f(DM , 
i=l 

where D x (i) refers to the i th generated sample. Simple Monte Carlo has two extremely useful properties. First, MC is 
guaranteed to be unbiased and thus, on average, will return the correct answer. Second, MC is guaranteed to converge 
to the correct answer at a rate of 0(n 1 ^ 2 ) independent of the number of dimensions . That is, for any problem, increasing 
the number of samples by a factor of one hundred decreases the expected error by a factor of ten. As Simple Monte 
Carlo has zero bias, all of the error comes from the variance of Monte Carlo runs. This variance is due a different kind 
of variance concerning f(x) itself; fluctuations in the Monte Carlo estimate are directly related to fluctuations in f(x). 
As a result, the smaller the variance of f(x), the smaller the expected error of Monte Carlo. 

Many different sample generation techniques exist to help reduce the variance of Monte Carlo. Importance sam- 
pling [1J generates samples from an alternate distribution q(x) (instead of the true distribution p(x)) and estimates 
integral ([T| as a weighted average of sample values 

/= J f(x)p(x)dx = J ^ X ^ X \ (x)dx , 

f(D,(i))p(DAi)) n . 
fis ~N^ q(D x (i)) ■ W 

Importance sampling is often used when the tails of a distribution have a measurable effect on /, but occur very 
infrequently. Using Simple MC, several million samples are needed to accurately capture a once-in-a-million event, 
but by using importance sampling unlikely outcomes can be made to occur more frequently reducing the total number 
of samples needed. 

Quasi-Monte Carlo (QMC) techniques reduce variance by choosing sample locations more carefully. Sample 
points generated from simple Monte Carlo will inevitably cluster in some locations and leave other locations void 
of samples. QMC methods are usually deterministic and reduce variance by spreading points evenly throughout the 
sample space. Two such methods are the scrambled Halton sequence |2| and the Sobol sequence 0. While often 
effective, due to deterministic sample generation they are not guaranteed to be unbiased. It is also difficult to generate 
points from an arbitrary p(x); most QMC algorithms generate sample points from a uniform distribution over the unit 
hypercube. 

2.3 Supervised Learning 

A second class of techniques for estimating integrals from data seeks to use the data samples more efficiently through 
the use of a fitting algorithm and supervised learning techniques. The fitting algorithm uses data samples to make 
a "fit" to the data, and the integral estimate is the integral of the fit. Examples of fitting algorithms include splines, 
high-order polynomials and Fourier expansions; even Simpson's rule computes the quadrature by fitting a piecewise 
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quadratic polynomial to the data samples. Fits incorporate the spatial distribution of sample points and thus often give 
more accurate estimates of / than MC. However, using a fitting algorithm can induce bias and may or may not exhibit 
convergence to the correct answer as the number of points increases. Additionally, when the number of data points is 
"too small" many fitting algorithms exhibit higher variance than MC, and, as a result, using a fitting algorithm can be 
worse than not using one, especially with low numbers of sample points and in high-dimensional spaces. It is usually 
impossible to know a priori how many points is "too small" making it difficult to know when to use a fitting algorithm. 

Additionally, it is difficult to know if a fit to the data samples is an accurate representation of f(x) at values of x not 
in the data set. Many fitting algorithms exhibit a property known as "overfitting" where a fit to the data is very accurate 
at x locations that are among the data samples, but is very inaccurate at x locations not among the data samples. As 
a result, one cannot use the data samples themselves to both make a fit and to evaluate its accuracy. The standard 
technique for addressing this issue is known as cross-validation, where a certain percentage of the data is used to make 
the fit and the rest of the data is used to evaluate its performance. Usually this process is repeated several times, and 
the best performing fit is used as the output of the fitting algorithm (a.k.a. winner-takes-all). 

Stacking is a technique first introduced by Wolpert [4] which improves upon cross-validation by combining the 
results of different fits. Wolpert describes stacking in his paper as " a strategy ... which is more sophisticated than 
winner-takes-all. Loosely speaking, this strategy is to combine the [fits] rather than choose one amongst them. ... 
[Stacking] can be viewed as a more flexible version of nonparametric statistics techniques like cross-validation, ... 
therefore it can be argued that for almost any generalization or classification problem ... one should use [stacking] 
rather than any single generalizer by itself." Interested readers can see stacking applied to improving regression 0, 
probability density estimation |6|, confidence interval estimation |7|, and optimization |8|. Stacking recently gained 
fame as being a major part of the first and second place algorithms in the Netflix competition |9|. 

In addition to combining the predictions of multiple fitting algorithms, stacking can also be used to improve the 
prediction of a single fitting algorithm. In this variant of stacking, one again partitions the full set of data D x into a 
subset S i used to make fits and the remainder of the data, and a subset S 2 = D x \ S 1 . However now one does not 
use S 2 to determine how to combine the predictions of multiple fitting algorithms that were all run on S 1 . Instead one 
uses S 2 to correct the prediction of a single algorithm that was run on S \. This variant of stacking is closest to the 
algorithm we use below. For a comparison of stacking to other generalizers, see Clarke[ 10|. 



3 Stacked Monte Carlo 

The main idea of Stacked Monte Carlo (StackMC) is to incorporate the variance reduction potential of a fitting al- 
gorithm while avoiding introducing bias and additional variance from poor fits. It makes no assumptions about the 
function f(x) nor the sample generation method, and therefore StackMC can be used to augment nearly any Monte 
Carlo application. 

Let us assume that we have a function g(x) that is a reasonable (though not necessarily perfect) fit to f(x). Equation 
([TJ can be re-written as 

/ = J ag(x)p(x) dx + J (f(x) - ag(x))p(x) dx 

ag + j (f(x) - ag(x))p(x) dx , (3) 

where a is a constant and g = J g(x)p(x) dx. Instead of using Monte Carlo samples to estimate ([T| directly, we can 
use the Monte Carlo samples to estimate the integral term in ([3]), i.e. 

1 N 

f^ag+-Yjf(D x (i))-ag(D x (i)) . (4) 

i=l 

Since g(x) is assumed to be a reasonable fit, for a properly chosen a, f(x) - ag(x) has lower variance than f(x) 
alone (see Fig. [TJ, and so a MC estimate of (|3]l has lower expected error than an estimate of ([T}. In fact, it can be 
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shown fin that the optimal value of a to minimize expected error is 

ff = p— , (5) 

where cry and cr g are the standard deviations of f(x) and g(x) respectively, and p is the correlation between / and g. 
Intuitively, if g is a good fit to /, p (and correspondingly a) will be high and g will be trusted as a good estimate for /. 
If g is a poor fit to /, p and a will both be low and g will be downplayed. Looking at it from a different perspective, 
Q takes the expected value estimated from g(x), and corrects it with a term estimating the bias of g(x). Either way, 
by using ([3]) the error should be lower than using either Monte Carlo or the fitting function alone. Since the expected 
value of the fit is both added and subtracted, (|4]) remains an unbiased estimate of fiji. Thus, Q incorporates a fitter 
while remaining unbiased, and a allows us to emphasize good fits while de-emphasizing poor ones. 

The obvious question is how to obtain g(x) and find a from data samples. The first step is to pick a functional 
form for a fitting algorithm g(x), such as a polynomial with unknown coefficients. g(x) can be nearly anything, the 
only restrictions are that it can make predictions at new x values. For now we also require that g can be analytically 
integrated over the volume, i.e. we can calculate 

I = J g(x)p(x) dx (6) 

analytically (see further discussion later in the paper). By comparing the output of a fit g(x) to the true f(x) values, 
an estimate for the "goodness" of the fit (and by extention a) could be obtained, and Q could be applied. However, 
doing this directly would cause over-fitting and would lead to an inaccurate estimate of a. 

Over- fitting can be mitigated by using a technique known as £-fold cross-validation IPO) . The N data samples 
are randomly partitioned into k testing sets, D x ' ,test ,i = l,...,k, which are mutually exclusive and contain the same 
number of samples (differing slightly if N/k is not an integer). Training sets, D x ,tram ,i = l,...,k, contain all of the 
data samples not in D x ' ,est . We call A^,- the number of samples in D x ,,est and m, the number of samples in D x l ' tram (such 
that Ni + ;«, = AO. One fit per fold, gi(x), is created using only the Ni samples in D x l,,ram (for a total of k fitters). The 
samples in a training set are called "held in" points because they are used to generate the fit, whereas points in the 
corresponding testing set are "held out" points because the fit is generated without using these samples. Each data 
sample is in a held out data set once and in the held in data set k — 1 times. Using gi(x) a prediction of the function 
values for the points in D x ' Jest is made ( gi{D x ' ,ram ) ). Standard cross-validation evaluates the accuracy of each of the 
fits, and chooses the best fit to use as a single g(x). 

Instead, we adopt the approach of stacking and use ([3) to get an estimate of / from each of the fitters, and we can 
use the held out data points as an estimate of the integral term. 
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fsmcii) =ag; + J ( f(x) - agi(x))p(x) dx (7) 
The mean of these estimates is taken as the final prediction 



1 - 

f ~ fsmc — T / .fsmci}) . (9) 

We can also use the predictions at the held out points to estimate a. crj is estimated from the variance of the 
sample function values themselves, cr g from the variance of the predictions, and p from the correlation between the 
predictions and the truth. 



1 N 



7=1 

k N, 



°7 



7=1 



ir iVi 



^ '=1 7=1 

1 * N ' 



i=l 7=1 



cov(/, g) 



Finally, we plug into (|9) 



i=i i=i ' j=i 



(10) 



One final correction is needed to complete StackMC. It is possible that some or all of the g differ greatly from / 
but a is calculated high because of good predictions at held out data points. If left alone, this would cause StackMC 
to return a low-quality prediction on some data sets. 

However, the error in the mean (EIM) of the Monte Carlo samples can be used as a second metric to evaluate the 
goodness of the fit. EIM is defined as 



\ Y 



and represents the uncertainty in the mean of a set of Monte Carlo samples. Specifically, it gives us a likelihood bound 
on / based on f mc and cry. We can find a "likelihood" for each fold by calculating 



Hi - fm 



(11) 
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The higher L,, the less likely it is that g ; = /, and the more likely it is that the fitter is bad and should be ignored. A 
heuristic is set that if max(L,) > C, f smc = f mc , i.e. all of the fits are ignored entirely and the Monte Carlo average is 
used. If C is set too low, then too many reasonable fits are ignored, and if C is set too high, then too many unreasonable 
fits are kept. A value of C = 5 was chosen based on experimentation; it was clear that setting C as low as 3 or as high 
as 7 were inferior. 

Finally, a note about the bias of StackMC. Because the expected value of the fit is being both added and subtracted, 
it is true that ^ (and by extension, |7]i and |9|) are also unbiased for a fixed a. However, using the held-out data to 
both set a and estimate the integral can introduce bias. In practice, we have found that this is only a problem for very 
small numbers of sample points where the output fit of the fitting algorithm changes significantly depending on the 
held-in samples. 



3.1 Generalized Stacked Monte Carlo 

The discussion above was for the case where the samples are generated according to a known p(x), this is only a 
special case of a class of estimation scenarios. While the overall methodology remains approximately the same for 
these other scenarios, some specifics (such as the calculation of a) change with the choice for g(x) and the sample 
generation method. 



3.1.1 Importance sampling 

If importance sampling methods are used, samples are generated from q(x) instead of p(x). As a result, ([T| is expanded 

as 

f = j ag{x)q(x)dx + 

and so g(x) should be a fitting algorithm for f ( - x)p( -^ instead of just f(x). As a result a few modifications are needed 
including the fact that, 

gi = J g(x)q(x) dx , 
and that in the calculation of a and ( |10) , — is used instead of just /. 



3.1.2 Unable to integrate g(x) over p(x) 

If the samples are generated from p(x), but the choice for g(x) cannot be integrated over p(x), there are two options. 
The first option is to use another function r(x) as a fitting algorithm for p(x) over which g(x) is integrable. We expand 
([T) as 



/f g(x)r(x) 
ag(x)r(x) dx + I fix) — a p(x) dx 
J ' P(x) 



Similarly to the above case, when calculating a and ( |10) , g is replaced with y. 

Alternately, when g(x) is significantly less computationally expensive than f(x), g itself can be estimated by Monte 
Carlo techniques. When estimating g from N g additional samples, it should be the case that N% » so that errors in 
the estimate of g do not cause significant errors in the estimate of /. 



3.2 A Simple Example 

We attempt to find the expected value of 3x 6 + 3.6x 5 -91.29x 4 - 19.41x 3 + 57.15x 2 - 14.43.x + 0.9, where x is generated 
according to a uniform distribution between -3 and 3. Thus, the integral of concern is 
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Left-out Fold 


X 


/(*) 


A> 




Pi 




gi(x) 


8i 


1 


0.4087 
-0.6950 
-0.0943 

0.1152 


0.2438 
7.8350 
0.9259 
-0.0166 


2.7385 


-4.3737 


-3.9500 


-9.4712 


-0.3549 
7.0498 
3.1237 
2.1675 


1.4281 


2 


0.4117 
0.2745 
0.1823 
0.2882 


0.2420 
0.1108 
-0.0163 
0.1342 


2.4683 


-3.3829 


-3.0183 


-11.3595 


-0.2284 
1.0774 
1.6825 
0.9708 


1 .4622 


3 


-0.6318 
-0.3923 
-0.8345 
0.7716 


7.6689 
4.7811 
6.4358 
-5.2874 


0.7900 


0.3755 


3.3397 


-22 0257 


7.4398 
2.4864 
15.6002 
-7.0489 


1.9032 


4 


0.5711 
-0.5988 

0.9607 
-0.6411 


-0.5683 
7.4412 
-16.6302 
7.7172 


1.8054 


-6.7386 


-0.2738 


-1.3468 


-2.3831 
6.0312 

-6.1154 
6.3677 


1.7141 


5 


0.7124 
0.1206 
-0.3960 
0.2816 


-3.2834 
-0.0208 
4.8377 
0.1230 


2.3965 


-3.8653 


-3.4306 


-10.9995 


-6.0740 
1.8609 
4.0722 
0.7901 


1.2529 



Table 1: Details of simple example calculations. 



i r 

f=2j (^ 6 + 1-2x 5 - 30.43jc 4 -6.47x 3 + 19.05x 2 -4.81x + 0.3)c/jc . 

The actual value of / can be calculated to be 0.7069. Twenty samples are generated and divided into five folds, 
with each fold's testing set containing four points and each corresponding training set containing the other sixteen 
points. While the function of interest is a sixth order polynomial, g(x) was chosen to be a third-order polynomial: 
Pq + [i\x + fax 2 + fox 3 . gi(x) is found by choosing the yS's which minimize the squared error of the data in D x l ' ,ram 
(using the standard pseudo-inverse). g2(x) - gs(x) are set similarly. Next, gi(D x '' tes '(j)) is evaluated for each test point 
in 

The following parameters are calculated to find a: 

(if = 1.1337 
H g = 1.9258 
07 = 5.6835 
a g = 5.2678 
cov(f,g) = 24.0999 
= covif^ = Qm9 

a =p— =0.8685 . 
Finally, equation ( fT0) i is applied to calculate f smc . 

Using the Monte Carlo estimate alone gives f mc = 1.1337, and using a fit to all of the data points alone gives 
ffit = 1.3412. By using StackMC, we find f smc = 0.8081, much closer to the true value than either of the individual 
estimates. Details of the calculations can be seen in Tables[T]and|2] and the fit from the first fold can be seen in Fig. [2] 
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Fold 




f(P*"U)) - a gi (D >■"■»(])) 


Corrected g ,• 


1 


1.4281 


-0.3553 


0.8795 


2 


1.4622 


-0.6428 


0.6271 


3 


1.9032 


-0.6122 


1.0407 


4 


1.7141 


-1.3569 


0.1318 


5 


1.2529 


-0.2732 


1.3613 


fsmc 


0.8081 



Table 2: Details of fold combination calculations. 




Figure 2: Example of the first fold. 
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4 Results 



We test the performance of StackMC on a number of different problems. In the first set of example cases the problems 
have known analytic results and an exact analysis of the performance of StackMC is examined. In the second set, 
StackMC is applied to two aerospace problems in the literature. While an exact solution to the aerospace problems 
is not available, an approximate answer is obtained from a Monte Carlo estimate using 100,000 samples. For each 
example problem, a range of number of samples was tested using two thousand simulations at each value of N. In 
each case, 10-fold cross-validation was used. Unless otherwise noted, the plots in this section have three lines which 
show the mean squared error versus the number of sample points. The lines represent the average squared error from 
using just the prediction of Monte Carlo (green), the fit to all the samples (red), and StackMC (blue). 

4.1 Analytic Test Cases 

4.1.1 Ten-dimensional Rosenbrock function under a uniform p(x), polynomial fitter. 

The D-dimensional Rosenbrock function is given by 

D-l 

i=i 

and x is drawn from a uniform distribution over the [-3,3] hypercube. 

The fitting algorithm is chosen as a third-order polynomial in each dimension whose form is 

D 

g(x) = /3 + ^PijXi +j3 2 ,iX? +Pxix] , 
;=i 

where the are free parameters that are set from the data samples. 

A comparison of the error of StackMC can be seen in Fig. [3] for the ten-dimensional version of the Rosenbrock 
function. At low numbers of sample points, Monte Carlo is more accurate, on average, than the polynomial fit to all 
the data points, but at higher numbers of points the polynomial outperforms Monte Carlo. Throughout the entire range 
of number of samples, StackMC performs at least as well as the best of the two. Additionally, as can be seen in Fig. [4] 
the polynomial fitter is actually a biased fitter for this example problem. StackMC is able to use cross-validation to 
remove the bias of the fitter while keeping the accuracy of its estimate. 

4.1.2 Ten dimensional Rosenbrock function under a Gaussian p(x). 

This is the same set-up as above, except that x is generated according to a multivariate Gaussian distribution, each 
dimension independent with p = and <x = 2. 

Much like the case above, in Fig. [5]it can be seen that at low numbers of sample points Monte Carlo outperforms 
the polynomial fit, whereas at high sample points the polynomial does better than Monte Carlo alone. However, for 
all numbers of sample points, StackMC outperforms the other two algorithms. 

4.1.3 20 Dimensional BTButterfly, uniform p(x), Fourier fitter 

In the previous examples, the Rosenbrock function is a A' 1 ' order polynomial, and the fitting algorithm is a 3 order 
polynomial, so StackMC had reasonable fits to use to improve upon Monte Carlo. In order to challenge StackMC, a 
function we call the BTButterfly was created so that it would be very difficult to fit accurately. 
Like the Rosenbrock function, its general form is given by 

D-l 

fix) = 2^h{xi,Xi+{) 

!=1 

with the contour plot of h shown in Fig. [6] x is generated uniformly over the [-3,3] hypercube. 
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10-D Rosenbrock Test Function with Uniform p(x) 
Average Squared Error vs. Number of Samples 
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Figure 3: Comparison of MC, fitter, and StackMC for the Rosenbrock function with uniform uncertainty. 
A Fourier expansion for g(x) was chosen whose form is given by 



g(x) =j3 + 2^/3ijCOs(Xi) + /3 2 j cos(2xi) + p Xi cos(3x,) +f3 4 jsm(x) i + /3 5 j sin(2x,) + f} 34 sin(3x,-) . 



!=1 

Results for this function are shown in FigJT| and exhibit the same trends seen previously; StackMC has as low of an 
error as the lowest of the two methods individually. 

4.2 Aerospace Applications 

4.2.1 Future Aircraft Uncertainty Quantification (UQ) 

In Ref.| 13 1, the authors use the Program for Aircraft Synthesis Studies [14| (PASS), a conceptual aircraft design tool, 
to predict the fuel burn of future aircraft given certain assumptions about technology advancement in the 2020 and 
2030 time frames. In their predictions for single-aisle aircraft in 2020, the authors model eight probabilistic variables 
representing different effects of improvements in aircraft technology (propulsion, structures, aerodynamics). Each of 
the variables is represented by a unique beta distribution. The authors generated 15,000 samples (each representing 
one optimized aircraft) from which they measured the expected fuel-burn metric and the standard deviation of the 
expected fuel-burn metric. 

To apply StackMC, we choose a third-order polynomial fitter. A polynomial fitter is convenient for this case 
because the E[x"] over a beta distribution is easily found analytically as follows: 



E[x n ] 



n? =1 <r+£ + i-l 



(12) 
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Figure 4: Expected output of MC, fitter, and StackMC with error in the mean for the 10-D Rosenbrock function. Note 
that the fitting function is biased, especially at lower numbers of sample points, while StackMC and MC are not. 
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10-D Rosenbrock Test Function with Gaussian p(x) 
Average Squared Error vs Number of Samples 

10 8 , , , , ■ ■ ■ 




Number of Samples 



Figure 5: Comparison of MC, fitter, and StackMC for the 10-D Rosenbrock function with Gaussian uncertainty. 
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Figure 6: Contour plot of successive dimensions of the BTButterfiy function. 
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20-D BTButterfly Test Function with uniform p(x) 
Average Squared Error vs Number of Samples 




Number of Samples 



Figure 7: Comparison of MC, fitter, and StackMC for the BTButterfly function with uniform uncertainty. 



where a and p are the two parameters of the beta distribution, B(a,/3) (not to be confused with the StackMC parameter 
a). 

A set of 100,000 samples were generated and the mean and standard deviation of their function values were taken as 
the "true" expected value and standard deviation. The standard deviation is defined as -\/E[x 2 ] - E[^] 2 , and therefore, 
to find the standard deviation we can run StackMC twice, once to find E[x] and a second time to find E[x 2 ]. The results 
of applying StackMC to find the expected value and standard deviation are shown in Fig.[8]and Fig. [9} respectively. 

The exact same trend is seen here as in all of the analytic problems. At low numbers of samples, where the fit to 
all the data points performs badly, StackMC does as well as Monte Carlo. At high numbers of samples, where the 
fit greatly outperforms Monte Carlo, StackMC does as well as the fitter. In between these two extremes, StackMC 
outperforms both algorithms. 

At 1500 samples, StackMC reduces the error by roughly an order of magnitude compared with the Monte Carlo 
result used by the authors. Each sample takes about 6 seconds to generate, whereas StackMC takes less than a half a 
second to form the 10 fits, calculate a, and calculate (JTOj . For a computational cost of less than one additional sample, 
StackMC has reduced the number of samples needed by 90% for the same expected error. 



4.2.2 Sonic Boom Uncertainty Quantification 

The uncertainty quantification of sonic boom noise is used as the second application case. Unlike the aircraft design 
test case, the response surface for the sonic boom noise signature is not smooth; the output can vary significantly with 
slight adjustments to the input parameters. Colonno and Alonso[15| recently created a new sonic boom propagation 
tool, SUBoom, and used it to analyze the robustness of several aircraft pressure signatures optimized for minimal 
boom noise. Specifically, in their high-fidelity near-field case, they have 62 uncertain variables: 4 representing aircraft 
parameters such as cruise Mach number and roll angle, and 58 representing uncertainties in the near-field pressure 
signal. Like in the Aircraft UQ case, 100,000 samples were generated and used as the true mean. A third-order 



polynomial was used as g(x). The results can be seen in Fig. 10 
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Aircraft UQ Test Case 
Average Squared Error vs Number of Samples 




Number of Samples 



Figure 8: Comparison of MC, fitter, and StackMC for the Aircraft UQ test case finding E[x]. 



Aircraft UQ E[x 2 ] 
Average Squared Error vs Number of Samples 




Number of Samples 



Figure 9: Comparison of MC, fitter, and StackMC for the Aircraft UQ test case finding E[x 2 ]. 



SUBoom Test Case 
Average Squared Error vs Number of Samples 



MC Estimate 




Number of Samples 

Figure 10: Comparison of MC, fitter, and StackMC for the SUBoom UQ test case. 

Even with the discontinuities in the function space, the same performance for StackMC is seen. StackMC does as 
well as the best of Monte Carlo and the fitting algorithm. The reduction in error is not as great here as in the Aircraft 
UQ case because a polynomial is not a good fit to f(x), though using domain knowledge to improve the accuracy of 
the fitting algorithm would further increase the performance of StackMC. Despite the relatively poor performance of 
the fitting algorithm, it is still better to use StackMC than to not regardless of the number of sample points used. In 
some senses, this test case demonstrates one of the strengths of StackMC: at virtually no additional computation cost 
the user is guaranteed the highest accuracy without having to worry about the appropriateness of the methodology for 
a given number of samples. 

5 Conclusion 

In this paper, we have introduced a new technique, Stacked Monte Carlo, which reduces error of Monte Carlo sampling 
by blending several different fitters of /(x). StackMC is unbiased, (thus retaining a major advantage of Monte Carlo), 
and is able to use cross validation to effectively incorporate the use of a fitting function. 

StackMC is an extremely generic post-processing technique. It is applied after the samples are generated and makes 
no assumptions about the generation method, and therefore StackMC can be used not only with simple MC, but also 
with other sample generation techniques such as Importance Sampling, Quasi-Monte Carlo, and Markov-Chain Monte 
Carlo. Furthermore, StackMC makes no assumptions about f(x); it not only applies to smooth functions, but also to 
discontinuous functions, and even functions with discrete variables. In a future paper, |[T6l we will present results 
of the application of StackMC to different input regimes and different sample generation methods. We will examine 
higher dimensional spaces, explore the application of StackMC to multi-fidelity methods, and extended StackMC to 
incorporate multiple fitting functions. 

The computation time of StackMC is dominated by forming the fits gi(x). As a result, StackMC is only affected 
by the curse of dimensionality insofar as the fitting algorithm is. In both of the aerospace applications, the additional 
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cost of the entire StackMC algorithm was virtually non-existant; there are significant improvements in accuracy for 
less computation time than the cost of one additional function evaluation. The only assumptions made about g(x) are 
that it be able to predict the value at new sample locations and we are able to evaluate g analytically (and even this 
second assumption is relaxed). As shown in the BTButterfiy and sonic boom example cases, the fit does not need to be 
particularly good to see improvement, and so the choice for g(x) can be modified as computational effort and domain 
knowledge allows. StackMC does not eliminate the need for finding better fitters; a better fitter will always lead to 
an improved result. With the exception of the error in the mean test (whose value we did not change for any of the 
example cases), StackMC requires no user-set parameters that need to be heuristically tuned to see good results. 

The version of StackMC implemented in this paper is relatively naive. Instead of taking the mean of the results 
from the different fits, taking a weighted combination of the fits could improve the estimate. StackMC currently only 
partitions the data into folds once, but we could imagine re-partitioning the same samples several times to have more 
fitters. We use Monte Carlo to estimate the integral term in (FT}, but using a fitter, or even using StackMC recursively, 
could improve the estimates of f smc . Furthermore, a was set as a constant, but in general a could vary between the 
folds or could even be a function of x so the confidence in the fit varies spatially. 

Despite this simplistic implementation, StackMC performs at least as well as the better of MC and the fitting 
function, and for a range of sample points outperforms both. Normally, there is a transition number of samples at 
which using a fit to all the data samples has a smaller average error than MC, and it is hard to know a priori if it is 
better to use the fit. StackMC eliminates the need to decide whether or not to use a fit; it will never be harmful to do 
so. While it is not true that for any set of data samples f smc is closer to / than f mc , on average StackMC reduces the 
expected error. StackMC is a very generic method for the post-processing of data samples; it can be used by anyone 
trying to estimate an integral or expected value of a function where p(x) is known. 
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